Parameter Estimation in Searches for the Stochastic Gravitational- Wave Background 
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The stochastic gravitational-wave background (SGWB) is expected to arise from the superposition 
of many independent and unresolved gravitational-wave signals of either cosmological or astrophys- 
ical origin. The spectral content of the SGWB carries signatures of the physics that generated it. 
We present a Bayesian framework for estimating the parameters associated with different SGWB 
models using data from gravitational-wave detectors. We apply this technique to recent results 
from LIGO to produce the first simultaneous 95% confidence level limits on multiple parameters in 
generic power-law SGWB models and in SGWB models of compact binary coalescences. We also 
estimate the sensitivity of the upcoming second-generation detectors such as Advanced LIGO/ Virgo 
to these models and demonstrate how SGWB measurements can be combined and compared with 
observations of individual compact binary coalescences in order to build confidence in the origin of 
an observed SGWB signal. In doing so, we demonstrate a novel means of differentiating between 
different sources of the SGWB. 

PACS numbers: 95.85.Sz, 97.60.Jd, 04.25. dg, 98.80.Cq 



Introduction. — The stochastic gravitational-wave 
background (SGWB) is expected to arise from the 
superposition of gravitational waves (GWs) from many 
uncorrelated and unresolved sources. Numerous cosmo- 
logical and astrophysical models have been proposed. 
Cosmological models include the slow-roll inflation 
model 0, 0], a variety of inflationary models with 
significant boosts in energy density at high frequencies 
(parametric resonance in the preheating phase Q, 
models of axionic natural inflation Q), models based 
on cusps or kinks in cosmic (super)strings [Hl-ll(ij|. and 
models of alternative cosmologies such as pre-Big-Bang 
models [ill 12 1- Astrophysical models integrate contri- 



butions from various astrophysical objects across the 
universe including compact binary coalescences (CBC) 
of binar y n eutron stars (BNS) or binary black holes 
fBBH) |l3-17j). magnetars 11814201 . rotating neutron 
stars I2l [-[25l . the first stars [26[, and white dwarf 
binaries [271 ]. 

In all models, the SGWB is described in terms of the 
normalized GW energy density, 

f_ dpcw 
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where rfpcw is the energy density of GWs in the fre- 
quency range / to / + df and p c is the critical energy 
density of the universe |28| . The amplitude and the fre- 
quency dependence of this GW spectrum depend on the 
physics of the model that generated it. For example, in 
the CBC model the spectrum is determined by the rate 
of binary systems throughout the universe and by the 
distribution of their chirp masses. 



Past SGWB searches, performed using data from the 
first-generation interferometric GW detectors LIGO 0, 
l30l | and Virgo [3l|, assumed a power-law model, 



^Gw(/)=^(///rof) Q , 



(2) 



and set limits only on the amplitude A for fixed values 
of the spectral index a and of the reference frequency 
/ref S3- 3f| . In this Letter we introduce a Bayesian tech- 
nique, building on [36| , to simultaneously estimate multi- 
ple parameters of SGWB models using cross-correlation 
data from pairs of GW detectors (H3] • We apply this tech- 
nique to recent results from LIGO to produce the first si- 
multaneous limits on multiple parameters for power-law 
and CBC models of the SGWB. We also estimate the 
sensitivity of second-generation GW detectors to these 
models — Advanced LIGO (aLIGO) 1371 1 . Advanced Virgo 
[H|, GEO-HF IH, and KAGRA [Sjllll are expected 
to produce first results in 2014, and will be sufficiently 
sensitive to probe a variety of interesting SGWB mod- 



els (see, e.g., [17j). We demonstrate that the technique 



can be easily extended to include other measurements, 
such as direct GW observations of resolvable compact 
binary coalescences, in order to better constrain astro- 
physical parameters and to gain insight into the origin of 
the observed SGWB signal. 

Method. Wc define the cross-correlation estimator 
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where T is the measurement time, s\(f) and S2(/) are 
Fourier transforms of the strain time-series of two GW 
detectors, and Q(f) is a filter J28[. In the small-signal ap- 
proximation, and assuming stationary, Gaussian noise — 
uncorrelated between the two detectors — the variance of 



2 



Y is given by: 



dfPi(f)P 2 (f)\Q(f)\' 



(4) 



where P%(f) are the one-sided strain power spectral den- 
sities of the two GW detectors. Optimization of the 
signal-to-noise ratio (SNR) leads to the following op- 
timal filter for a frequency-independent GW spectrum 
^gw(/) = HI]: 



7(f) 



(5) 



where 7(/) is the overlap reduction function arising from 
the overlap of antenna patterns of GW detectors at dif- 
ferent locations and with different orientations [28| . The 
normalization constant M is chosen so that (Y) = fio- 

Most proposed SGWB models flu{f\9) vary slowly 
with frequency; (here 9 denotes the model parameters). 
It is therefore a good approximation to compute the 
above estimator for a series of small, 0.25Hz-wide fre- 
quency bins (see, e.g., (HI), and define the following like- 
lihood function: 



L(Yi, <jj|0) oc exp 



9 2^ rfi 



(6) 



Yi and <Ji 34[ to place 95% confidence level (CL) limits 
on a and A. These limits are the first to simultaneously 
constrain two parameters of an SGWB model. We also 
calculate the projected sensitivity for colocated aLIGO 
detectors that would be obtained after 1 yr of exposure at 
the design strain sensitivity and with SNR = 2 (assuming 
Y t = 0). 

Compact binary coalescences. — Coalescences of binary 
systems such as BNS and BBH are among the most 
promising sources of GWs in the 10 — 1000 Hz frequency 
band. The nearest CBCs are expected to produce GW 
signals strong enough to be individually detected by 
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where the sum runs over frequency bins fi, and Yi and 
cf are the estimator and variance in the frequency bin 
fi given respectively by Eqs. |3] and 01 Since Y is cal- 
culated by averaging over many short time segments, its 
distribution is expected to be Gaussian due to the cen- 
tral limit theorem, as observed, e.g., in j34|. Multiplying 
the likelihood with the prior distribution for 9 yields the 
Bayesian posterior distribution for the free parameters 9, 
which can then be used to extract confidence intervals 
for 9. In the subsequent results, we take all priors to be 
flat within the plotted range of 9 (and zero elsewhere) . 

Power-law spectrum. — As most SGWB models predict 
power-law dependence in the LIGO frequency band, we 
apply the above formalism to the power-law model of the 
SGWB (see Eq. [5J with a reference frequency of / ro f = 
100 Hz. In Fig. [1] we use recent LIGO measurements of 



FIG. 1: The solid curve denotes the 95% CL limit for the 
power-law SGWB model Q M (f;A,a) = A(//100 Hz) Q , in the 
a- A plane using the latest SGWB measurement with LIGO 
detectors [13]. The dashed curve denotes the projected sensi- 
tivity obtainable by aLIGO at SNR = 2. assuming 1 yr of data 
from two colocated detectors operating at design sensitivity 
(assuming Yi =0). 



second-generation detectors [42j. Meanwhile, integrat- 
ing contributions from all CBCs in the universe leads to 
an SGWB that may also be detected by cross-correlating 



pairs of GW detectors 



Following [17j , in the BNS case we include only the in- 
spiral part of the GW signal, and write the GW spectrum 
as an integral over redshift z and over the chirp mass M' c : 
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We have verified that the shape of the chirp mass dis- 



tribution has a negligible effect on the GW spectrum, 



3 



especially at frequencies below the peak where terrestrial 
GW detectors are most sensitive (see Fig. [2|. We there- 
fore approximate the GW spectrum using only the aver- 
age chirp mass M c , as shown in the second line of Eq. 
[7J The free parameters of the model are therefore M c 
and A, the mass fraction parameter. We treat the BNS 
and BBH populations separately, with different average 
chirp masses. Furthermore, for the BBH case we use the 
more complex functional form derived by and used 
by Ifljllll, which includes the inspiral, merger, and ring- 
down contributions to the gravitational-wave signal (see 
[lij for more detail). The merger and ringdown spectra 
are computed assuming equal mass black holes. 

The A parameter is proportional to the local CBC rate 
per unit volume (discussed further below), and it cap- 
tures uncertainties associated with the mass fraction of 
the neutron star and black hole progenitors, the fraction 
of massive binaries formed among all stars, and the frac- 
tion of all binaries that remain bounded after the second 
supernova event. Further, G is Newton's constant, c is 
the speed of light, E(n m , Q\, z) = ^/fi m (l + z) 3 + 
captures the dependence of the comoving volume on red- 
shift (we use the standard ACDM cosmology, fl m = 0.3, 
ft A = 0.7), and Rv{z) is the observed rate of binary co- 
alescences given by the following integral over the time- 
delay td between creation and coalescence of the binary: 



R v (z) 



1 



z f 



R*(t c (z)-t d )p(t d )dt d . (8) 



Here t c (z) is the cosmic time to coalescence correspond- 
ing to redshift z, Zf is the redshift at the formation time 
t c (z) — td, and i?* is the star formation rate (SFR). Mul- 
tiple SFR models have been proposed |45|-|49[ . We adopt 
the model from [ifl and note that other SFR models can 
lead to variations in J7m(/; M c , A) of up to a factor of two 



17j. Finally, p(td) is the probability density function for 
the time-delay — we use p(t d ) ~ t~{ , with the minimum 
time-delay i min = 20 Myr for BNS and 100 Myr for BBH 
cases, and with f max equal to the age of the universe. We 
note that other choices of time-delay distribution could 
also lead to a factor of two difference [l7|. Finally, note 
the local CBC rate is given by i?i oc = Ai?y(0). 

The upper limit on the integral range in Eq. [JJ depends 
on both the emission frequency range, / m ; n — / max , in the 
source frame, and on the maximum redshift z max = 6 
considered for the star formation history calculation: 



^max 



//-i 



if / < /max/(l + Zmax) 

otherwise 



(9) 



The top-left panel of Fig. [5] shows several examples of 
BNS and BBH spectra — depending on the choice of pa- 
rameter values, the spectra could peak in the sensitive 
band of aLIGO. In the top-right panel of Fig. [2] we show 
the 95% CL limits in the i?i oc — M c plane (or, equiv- 
alentl y, A — M c plane) obtained using recent LIGO re- 
sults [34| and the projected sensitivities for aLIGO at 



SNR = 2, assuming lyr of exposure for colocated de- 
tectors operating at design sensitivity (and % = 0) [HBj ]. 
Note that second-generation detectors are expected to be 
10 4 x more sensitive to the CBC SGWB models than ini- 
tial LIGO [34j. Also note that, as stated above, different 
choices of the star formation rate and time-delay distri- 
bution lead to a factor of 2 uncertainty in the computed 
upper limits and projected sensitivities. 

The bottom row of Fig. [2] shows projected confidence 
contours for plausible simulated BNS (left) and BBH 
(right) signals that are within the reach of advanced 
detectors. For the BBH case, we consider the local 
BBH rate R loc = 0.17Mpc~ 3 Myr _1 , which is in be- 
tween the realistic rate (0.005 Mpc _3 Myr _1 ) and the 
optimistic rate (0.3 Mpc _3 Myr _1 ) from 42j. This re- 
sults in relatively narrow contours in the i?i oc — M c 
plane. It is evident, however, that this signal is not 
strong enough to break the degeneracy between the i?i oc 
and M c parameters — roughly 10 times higher rate of 
BBH systems is needed in order to break this degen- 
eracy. For the BNS case, we consider the local BNS 
rate i?i oc = 2 Mpc _3 Myr _1 , which is in between the 
realistic rate (1 Mpc~ 3 Myr" 1 ) and the optimistic rate 
(10 Mpc _3 Myr _1 ) from [42j. Again, the stochastic mea- 
surement alone is not sufficient to break the degener- 
acy between i?i oc and M c . However, we demonstrate 
in this case that our likelihood formalism can be nat- 
urally extended to include measurements of individually 
resolvable binary coalescences to extract more informa- 
tion about these parameters and about the source of the 
GW background. 

The R\ oc -M c likelihood function from individual CBC 
detections is given approximately by 



L(h, Al-Riocj M c ) cx exp 



(A - M c f 



(V(M c )R\ oc T) n V(M c )R loc T 

77.1 



(10) 



where fi is the average reconstructed chirp mass (with un- 
certainty cr^), h is the number of CBC detections above 
threshold, V(M C ) is the average volume of space in which 
CBC signals with chirp mass M c can be detected, and T 
is the observation time. Eq.rT0]does not take into account 
the fact that the detectability of CBC signals varies with 
chirp mass — SNR grows like M^ 6 [H(| — but this intro- 
duces only a small bias < 4%. The uncertainty de- 
pends on the variance of the true chirp mass distribution, 
the number of measurements, and the measurement er- 
rors associated with individual M c measurements, which 
are typically < 0.5M Q depending on the SNR. 

We assume that the true chirp mass distribution 
is a delta function at /i = 1.22M© (each NS has a 
mass of I.4M0). Given our "semi-optimistic" rate of 
2 Mpc _3 Myr _1 , we expect 80 BNS events in 1 year 
of aLIGO data 421. We further assume that each 
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chirp mass can be reconstructed with a precision of 
ps O.5M so that the uncertainty for the average is 
(7 M = 0.5MQ/V80 ps O.O6M . (If the true chirp mass 
distribution has a width comparable to or greater than 
the detector resolution then will grow appreciably.) 

In the bottom-left panel of Fig.[2]we include 95% con- 
fidence contours obtained using just stochastic data, just 
direct BNS observations, and the contour obtained by 
computing the joint likelihood function as the product 
of likelihood functions given in Eqs. [6] and QUI The 
joint likelihood leads to the tightest constraints on the 
model parameters, although in this particular example, 
the contour is determined primarily by the direct BNS 
observations. 

Perhaps an even more important application of the 
joint likelihood technique will be in understanding the en- 
ergy budget of the observed SGWB. In our example, the 
agreement between the stochastic and CBC contours in- 
dicates that the observed SGWB is consistent with a sin- 
gle SGWB source, namely the BNS coalescences. In gen- 
eral, however, multiple SGWB sources (BNS, BBH, mag- 
netars, and others) are expected to contribute to the ob- 
served SGWB. By constructing a "grand likelihood func- 
tion" with stochastic and transient measurements — and 
allowing for multiple SGWB models — it will be possible 
to identify contributions from different SGWB sources. 
Understanding the energy budget of the observed SGWB 
may be critical for detection of the cosmological contribu- 
tion, which will likely be masked by astrophysical SGWB 
sources in the frequency band of terrestrial GW detec- 
tors. Successful implementation of this technique will 
require simultaneous estimation of more than two pa- 
rameters, which in turn will require more sophisticated 
techniques for likelihood maximization, such as Markov- 
chain Monte Carlo 5l| or nested sampling [52|. These 
approaches are currently being investigated. 

Conclusions. — We have presented the first multi- 
dimensional confidence intervals for models of the 
stochastic gravitational- wave background. The two cases 
we considered, generic power-law and compact binary co- 
alescence models, can both be parametrized with just 
two parameters. Other models, invoking, for example, 
magnetars or cosmic strings, are likely to require more 
than two parameters. As the dimensionality of the likeli- 
hood function increases, techniques such as Markov-chain 
Monte Carlo 
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or nested sampling [52| will likely be 
necessary in order to efficiently calculate confidence in- 
tervals on model parameters. 

We have demonstrated that our technique can be ex- 
tended to combine stochastic measurements with direct 
detections of gravitational-wave transients, yielding bet- 
ter constraints on astrophysical parameters. Further- 
more, the technique provides a formalism for under- 
standing the energy budget of the observed stochas- 
tic gravitational-wave background, identifying contribu- 
tions due to different known astrophysical and cosmolog- 
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FIG. 2: Top-left: GW spectra for several examples of CBC models are shown. The solid and dot-dashed curves were computed 
using the model described here, for the BNS (black) and BBH (gray) cases, while the dashed curves were computed using the 
distribution of chirp masses from StarTrack simulation with the same average chirp masses as those used for the solid curves. 
The solid and dashed curves are nearly identical for low values of M c , indicating that using the average chirp mass is sufficient 
for the analysis presented here. Top-right: 95% CL limits for CBC models for initial LIGO [33] and projected sensitivities 
for aLIGO assuming SNR = 2 and 1 yr of coincident data from colocated interferometers operating at design sensitivity (and 
Yi = 0). The region M c < 2.5M Q corresponds to BNS models, and M c > 2.5M Q to BBH models. The three different styles 
(solid, dashed, and dot-dashed) correspond to three different choices of star formation rate R*(z) and time-delay probability 
density function p(td), indicating the uncertainty due to these choices at the level of a factor of 2 [l7(| . Bottom-left: projected 
95% CL contours for a simulated BNS signal measured with aLIGO. The simulated signal amplitude corresponds to local BNS 
rate in between the "realistic" and "optimistic" rates described in [Z^]. Combining the stochastic measurement with direct 
observations of individual BNSs leads to a marginally tighter constraint on the astrophysical parameters describing the BNS 
population. Bottom-right: projected 95% CL contours for a simulated stochastic BBH signal measured with aLIGO. A rather 
optimistic local rate R\ oc — 0.17 Mpc~ 3 Myr~ 1 was assumed for this simulation [4^], resulting in a relatively narrow contours 
in the R\ oc — M c plane. However, the signal is still not strong enough to differentiate between a high i?i oc - low M c signal and 
a low _Ri oc - high M c signal. 
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